library(tidyverse)
library(terra)
library(sf)
library(mapview)
library(raster)
library(rgeos)
library(stars)
library(lubridate)
library(ggplot2)
library(exactextractr)
library(patchwoRk)
library(gridExtra)USDA National Forest Type Group Dataset (2004): -Douglas-Fir, Ponderosa Pine, Fir-Spruce-Mountain Hemlock, Lodgepole Pine
# forest type groups and key
conus_forestgroup <- raster('data/conus_forestgroup.tif')
forest_codes <- read_csv('data/forestgroupcodes.csv')
# set crs
crs = crs(conus_forestgroup)EPA level-3 Ecoregions: -Canadian Rockies, Idaho Batholith, Middle Rockies, Southern Rockies, Columbian Mountains - Northern Rockies, Uinta Mountains
# level 3 ecoregions
l3eco <- st_read('data/us_eco_l3.shp') %>%
st_transform(., crs=crs)## Reading layer `us_eco_l3' from data source
## `G:\Other computers\My Laptop\Documents\Grad School\Research\WesternConiferRegeneration\data\us_eco_l3.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1250 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -2356069 ymin: 272048.5 xmax: 2258225 ymax: 3172577
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic_USGS_version# select northern rocky mountains from level3 ecoregions
eco_select <- l3eco %>%
filter(NA_L3NAME %in% c('Canadian Rockies','Columbia Mountains/Northern Rockies','Middle Rockies','Idaho Batholith'))
mapview(eco_select)MTBS Fires: -1988-1996 -5000+ acres of burning -500+ acres of high-severity -Within selected ecoregions ->25% of selected forest types
# mtbs fire perimeters
mtbs_full <- st_read('data/mtbs_perims_DD.shp') %>%
st_transform(., crs=crs)## Reading layer `mtbs_perims_DD' from data source
## `G:\Other computers\My Laptop\Documents\Grad School\Research\WesternConiferRegeneration\data\mtbs_perims_DD.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 29533 features and 22 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -166.1885 ymin: 17.94736 xmax: -65.33821 ymax: 70.15893
## Geodetic CRS: NAD83# select fires qualities of interest
# region of interest, 1988-1996, over 5k acres, some high severity
mtbs_select <- mtbs_full %>%
mutate(state = str_sub(Event_ID,0,2),
year = year(as.Date(Ig_Date))) %>%
filter(state %in% c("WA","ID","MT","WY","SD"),
Incid_Type == "Wildfire",
between(Ig_Date, as.Date('1988-01-1'), as.Date('1996-12-31'))) %>%
group_by(Event_ID,Incid_Name,Ig_Date,year) %>%
summarize(geometry= st_union(geometry),
acres=sum(BurnBndAc),
min_High_T = min(High_T),
max_High_T = max(High_T)) %>%
filter(min_High_T != -9999,
max_High_T != 9999,
acres >= 5000) ## `summarise()` has grouped output by 'Event_ID', 'Incid_Name', 'Ig_Date'. You
## can override using the `.groups` argument.# assign proportions of forest type to each fire polygon
fires_join <- st_join(mtbs_select,eco_select,join=st_intersects,left=FALSE) %>%
left_join(., exact_extract(conus_forestgroup,mtbs_select, append_cols = TRUE, max_cells_in_memory = 3e+08,
fun = function(value, coverage_fraction) {
data.frame(value = value,
frac = coverage_fraction / sum(coverage_fraction)) %>%
group_by(value) %>%
summarize(freq = sum(frac), .groups = 'drop') %>%
pivot_wider(names_from = 'value',
names_prefix = 'freq_',
values_from = 'freq')}) %>%
mutate(across(starts_with('freq'), replace_na, 0)))##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%## Joining, by = c("Event_ID", "Incid_Name", "Ig_Date", "year", "acres",
## "min_High_T", "max_High_T")# remove duplicates from ecoregion boundary crossover
# remove unnecessary columns, cleanup names
# filter to remove majority other forest types, majority not-forested cover
fires <- fires_join %>%
dplyr::select("Event_ID","Incid_Name","year","Ig_Date","acres","US_L3NAME","freq_0","freq_200","freq_220","freq_260","freq_280") %>%
rename("fire_name"="Incid_Name",
"fire_acres"="acres",
"date"="Ig_Date",
"ecoregion" = "US_L3NAME",
"freq_df"="freq_200",
"freq_pp"="freq_220",
"freq_fs"="freq_260",
"freq_lpp"="freq_280") %>%
mutate(freq_allother = 1-(freq_0 + freq_df+freq_pp+freq_fs+freq_lpp),
freq_forested = 1- freq_0,
freq_ideal = freq_df+freq_pp+freq_fs+freq_lpp)%>%
mutate(across(starts_with('freq'), round,2)) %>%
.[!duplicated(.[1]),]%>%
filter(freq_ideal >= 0.25)# import all mtbs rasters via a list
rastlist <- list.files(path = "data/MTBSmosaic", pattern='.tif', all.files=TRUE, full.names=TRUE)
allrasters <- lapply(rastlist, raster)
names(allrasters) <- str_c("y", str_sub(rastlist,28,31))
# create empty dataframe
severity_list <- list()
# loop through mtbs mosasics for 1986-1996
# extract burn severity raster for all selected fires
# calculate burn severity percentages for each fire
for (i in names(allrasters)){
mtbs_year <- allrasters[[i]]
fire_year <- filter(fires, year==str_sub(i,2,5))
raster_extract <- exact_extract(mtbs_year,fire_year, max_cells_in_memory = 3e+09,coverage_area=TRUE)
names(raster_extract) <- fire_year$Event_ID
output_select <- bind_rows(raster_extract, .id = "Event_ID")%>%
group_by(Event_ID , value) %>%
summarize(total_area = sum(coverage_area)) %>%
group_by(Event_ID) %>%
mutate(proportion = total_area/sum(total_area))%>%
dplyr::select("Event_ID","value","proportion") %>%
spread(.,key="value",value = "proportion")
severity_list[[i]] <- output_select
}## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:NA)##
|
| | 0%
|
|== | 3%
|
|==== | 5%
|
|====== | 8%
|
|======== | 11%
|
|========= | 14%
|
|=========== | 16%
|
|============= | 19%
|
|=============== | 22%
|
|================= | 24%
|
|=================== | 27%
|
|===================== | 30%
|
|======================= | 32%
|
|========================= | 35%
|
|========================== | 38%
|
|============================ | 41%
|
|============================== | 43%
|
|================================ | 46%
|
|================================== | 49%
|
|==================================== | 51%
|
|====================================== | 54%
|
|======================================== | 57%
|
|========================================== | 59%
|
|============================================ | 62%
|
|============================================= | 65%
|
|=============================================== | 68%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|===================================================== | 76%
|
|======================================================= | 78%
|
|========================================================= | 81%
|
|=========================================================== | 84%
|
|============================================================= | 86%
|
|============================================================== | 89%
|
|================================================================ | 92%
|
|================================================================== | 95%
|
|==================================================================== | 97%
|
|======================================================================| 100%## `summarise()` has grouped output by 'Event_ID'. You can override using the
## `.groups` argument.## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:NA)##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%## `summarise()` has grouped output by 'Event_ID'. You can override using the
## `.groups` argument.## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:NA)##
|
| | 0%
|
|=================================== | 50%
|
|======================================================================| 100%## `summarise()` has grouped output by 'Event_ID'. You can override using the
## `.groups` argument.## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:NA)##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%## `summarise()` has grouped output by 'Event_ID'. You can override using the
## `.groups` argument.## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:NA)##
|
| | 0%
|
|============ | 17%
|
|======================= | 33%
|
|=================================== | 50%
|
|=============================================== | 67%
|
|========================================================== | 83%
|
|======================================================================| 100%## `summarise()` has grouped output by 'Event_ID'. You can override using the
## `.groups` argument.## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:NA)##
|
| | 0%
|
|==== | 6%
|
|======== | 11%
|
|============ | 17%
|
|================ | 22%
|
|=================== | 28%
|
|======================= | 33%
|
|=========================== | 39%
|
|=============================== | 44%
|
|=================================== | 50%
|
|======================================= | 56%
|
|=========================================== | 61%
|
|=============================================== | 67%
|
|=================================================== | 72%
|
|====================================================== | 78%
|
|========================================================== | 83%
|
|============================================================== | 89%
|
|================================================================== | 94%
|
|======================================================================| 100%## `summarise()` has grouped output by 'Event_ID'. You can override using the
## `.groups` argument.## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:NA)##
|
| | 0%
|
|============== | 20%
|
|============================ | 40%
|
|========================================== | 60%
|
|======================================================== | 80%
|
|======================================================================| 100%## `summarise()` has grouped output by 'Event_ID'. You can override using the
## `.groups` argument.# combine extracted raster datasets
severity_df <- do.call(rbind, severity_list)
# join burn severity % to fires polygons
# fix naming
# filter dataset for 500 acres high severity
fires_severity <- left_join(fires,severity_df,by="Event_ID")%>%
rename(noburn= "1",lowsev = "2", medsev = "3", highsev = "4",regrowth = "5", error = "6") %>%
dplyr::select(- "NaN",-"regrowth",-"error") %>%
mutate(highsev_acres = fire_acres*highsev) %>%
filter(highsev_acres > 500) %>%
.[!duplicated(.[1]),]
# fix that one fire
# fires_severity[42,1] <- "FALLS 2"
# fires_severity[61,1]<- "HELEN CREEK 2"
# fires_severity[146,1]<- "TRAIL CREEK 2"
fires_select <- fires_severity %>%
left_join(.,exact_extract(conus_forestgroup,fires_severity, 'mode', append_cols = TRUE, max_cells_in_memory = 3e+08)) ##
|
| | 0%
|
|= | 1%
|
|== | 3%
|
|=== | 4%
|
|==== | 5%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============== | 20%
|
|=============== | 21%
|
|================ | 22%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|============================ | 39%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|=============================== | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|===================================== | 53%
|
|====================================== | 54%
|
|======================================= | 55%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================= | 70%
|
|================================================== | 71%
|
|=================================================== | 72%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|=============================================================== | 89%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 99%
|
|======================================================================| 100%## Joining, by = c("Event_ID", "fire_name", "year", "date", "fire_acres",
## "ecoregion", "freq_0", "freq_df", "freq_pp", "freq_fs", "freq_lpp",
## "freq_allother", "freq_forested", "freq_ideal", "noburn", "lowsev", "medsev",
## "highsev", "highsev_acres")fires_select$mode <- as.factor(fires_select$mode)
fires_select <- fires_select %>%
mutate(fire_foresttype = case_when(mode==200 ~ "Douglas-Fir",
mode==220 ~ "Ponderosa",
mode==260 ~ "Fir-Spruce",
mode==280 ~ "Lodegepole Pine",
TRUE ~ "Other"))
fires88 = fires_select %>% filter(year ==1988)Selected fires by year
# plot
mapview(fires_select, zcol = "year")Selected fires by majority forest type
# plot
mapview(fires_select, zcol = "fire_foresttype")